###equal-mass discretization of log normal
function discretize_lognormal(μ::Float64, σ::Float64, n::Int64)
    dist = LogNormal(μ, σ) #set up log normal distribution
    shocks = zeros(n) #preallocate market luck states

    for i = 1:n #begin loop to fill discretized vector
        quant = (2*i-1)/(2*n) #desired quantile
        shocks[i] = quantile(dist, quant) #fill element of state vector
    end
    shocks #return deliverable
end

###equal mass discretization of normal
function discretize_normal(μ_a::Float64, σ_a::Float64, n::Int64)
    dist = Normal(μ_a, σ_a) #set up log normal distribution
    draws= zeros(n) #preallocate market luck states

    for i = 1:n #begin loop to fill discretized vector
        quant = (2*i-1)/(2*n) #desired quantile
        draws[i] = quantile(dist, quant) #fill element of state vector
    end
    draws #return deliverable
end

###equal-mass distribution of joint normal distribution b/t parent HC child ability
function discretize_biv_normal_cond(μ_a::Float64, σ_a::Float64, μ_h::Float64, σ_h::Float64, ρ_ha::Float64, h::Float64, n::Int64)
    mean = μ_a + ρ_ha*(σ_a/σ_h)*(h - μ_h)
    var = (1-ρ_ha^2)*σ_a^2

    dist = Normal(mean, sqrt(var)) #set up log normal distribution
    draws= zeros(n) #preallocate market luck states

    for i = 1:n #begin loop to fill discretized vector
        quant = (2*i-1)/(2*n) #desired quantile
        draws[i] = quantile(dist, quant) #fill element of state vector
    end
    draws #return deliverable
end

###utility function
function utility(c::Float64)
    util=-1e10
    if c>=0
        util = log(c)
    end
    util
end

#function for finding index equivalent of a value in HC grid
function findmatch(index::Float64, target::Float64, interp::Interpolations.BSplineInterpolation)
    val = abs(interp(index)-target)
    val #return
end

#function for getting interpolated index of a given Float with an arbitrary grid, assuming linear interpolation
function get_index(val::Float64, grid::Array{Float64,1})
    n = length(grid)
    index = 0 #preallocation
    if val<=grid[1] #LEQ smallest element
        index = 1
    elseif val>=grid[n] #GEQ biggest element
        index = n
    else
        index_upper = findfirst(z->z>val, grid)
        index_lower = index_upper - 1
        val_upper, val_lower = grid[index_upper], grid[index_lower] #values
        index = index_lower + (val - val_lower)  / (val_upper - val_lower) #weighted average
    end
    index #return
end

#standard normal CDF
function 𝚽(x::Float64)
    return cdf(Normal(), x)
end

#mills ratio
function Mills(x::Float64, μ::Float64, σ::Float64)
    frac = (x - μ)/σ
    mills = μ + σ * pdf(Normal(), frac)/(1-cdf(Normal(), frac))
    mills #return
end


##grant grant calculator
function Grants(I::Float64, ap::Int64, sched::Array{Float64,2})
    grant, grant_fed=0, 0

    index_ability = 1
    if ap == 2 || ap == 3 || ap == 4
        index_ability = 2
    elseif ap==5
        index_ability = 3
    end

    #starting row
    row_start = 0
    row_select = row_start + 1
    for i = 2:14 #loop over rows
        if I>sched[row_start + i-1, 1] && I<sched[row_start + i, 1]
            row_select = row_start + i
        end
    end
    grant = sched[row_select, 1 + index_ability] #update grant based on ability
    grant_fed = sched[row_select, 5]
    grant_total = grant + grant_fed
    grant_total
end
